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EXECUTIVE  SUMMARY 


OBJECTIVES 

To  develop  a  mathematical  model  for  a  transducer  array  which  can  handle  nonlinear 
devices  such  as  the  HLF-6A  hydraulic  transducers.  The  model  should  be  capable  of 
handling  both  sinusoidal  and  nonsinusoidal  electrical  inputs. 


APPROACH 

Since  the  model  must  be  able  to  handle  nonlinear  devices  it  is  necessary  to  formulate 
the  model  in  the  time  domain  rather  than  in  the  frequency  domain.  The  internal  com¬ 
ponents  of  the  transducers  were  modeled  using  standard  hydraulic  flow  equations  which 
are  in  the  form  of  nonlinear  ordinary  differential  equations.  The  equations  are  nonlinear 
since  the  pressure  drop  across  a  hydraulic  valve  is  proportional  to  the  square  of  the  flow 
through  the  valve.  The  acoustic  radiation  interaction  impedances  were  obtained  in  the 
frequency  domain  using  the  boundary  element  program  CHIEF.  The  frequency-domain 
interaction  impedances  were  then  transformed  to  the  time  domain  using  an  FFT  routine 
to  produce  interaction  impulse  responses.  The  radiation  impedance  relations  in  the  fre¬ 
quency  domain  become  convolution  relations  in  the  time  domain.  When  these  radiation 
convolution  integral  relations  are  combined  with  the  model  for  the  transducer  internal 
components,  the  result  is  a  coupled  set  of  nonlinear  differential-integral  equations  that 
is  solved  numerically  by  an  explicit  time-stepping  procedure  similar  to  a  second-order 
Runge-Kutta  method. 

RESULTS 

Our  time-domain  model  was  first  applied  to  an  array  of  simple  linear  transducers  for 
which  the  solution  could  be  obtained  by  conventional  frequency-domain  methods.  The 
agreement  between  the  results  of  our  time-domain  model  and  the  frequency-domain  results 
was  excellent.  A  comparison  was  also  made  between  our  model  and  a  model  developed  by 
Hydroacoustics,  Inc.  for  a  single  transducer  element.  Again  the  agreement  was  quite  good. 
Finally,  predicted  results  using  our  time-domain  model  were  compared  with  experimental 
results  for  a  three-element  close-packed  array.  The  agreement  was  reasonably  good  in 
both  level  and  trends. 
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1  INTRODUCTION 


This  report  contains  a  description  of  a  mathematical  model  for  a  hydroacoustic  trans¬ 
ducer  array.  Although  the  model  was  developed  with  the  HLF-6A  transducer  in  mind,  it 
can  be  applied  to  other  hydroacoustic  transducer  arrays.  This  model  is  based  partly  on 
previous  modelling  work  for  the  HLF-4  transducer  [1],  but  involves  new  work  relating  to 
acoustic  radiation  loading.  Previous  models  have  attempted  to  couple  frequency-domain 
models  for  the  acoustic  field  with  time-domain  models  for  the  internal  components  of 
the  transducer.  In  this  report  we  present  a  complete  time-domain  model  for  the  coupled 
transducer/ acoustic-field  problem  which  is  applicable  to  both  single  elements  and  arrays. 
This  model  can  also  handle  nonsinusoidal  signals,  which  is  important  since  linear  superpo¬ 
sition  is  not  valid  for  nonlinear  devices.  The  acoustic  radiation  model  involves  obtaining 
certain  interaction  coefficients  in  the  frequency  dom^n  with  a  boundary  element  program 
called  CHIEF  [2,3]  and  then  transforming  these  coefficients  to  the  time  domain  by  using 
Fourier  transforms.  When  these  coefficients  are  combined  with  a  model  for  the  internal 
components,  the  result  is  a  coupled  set  of  differential-integral  equations  that  can  be  solved 
by  an  explicit  time-stepping  procedure  with  prescribed  initial  conditions. 

Since  the  model  for  acoustic  radiation  loauiing  is  new,  it  was  first  applied  to  a  simple 
linear  transducer  array  for  which  the  solution  could  also  be  obtained  by  the  standard 
frequency-domain  approach.  Section  2  contains  a  description  of  the  linear  model  that 
includes  the  new  acoustic  radiation  approach.  Results  of  the  time-domain  model  are  com¬ 
pared  with  corresponding  frequency-domain  results  in  order  to  validate  the  new  acoustic 
radiation  approach.  The  agreement  was  excellent.  Section  3  contains  a  description  of 
the  models  for  the  various  components  of  the  HLF-6A  hydraulic  transducer,  such  as  the 
first-stage  servovalve,  the  second-stage  valve,  the  drive  cavity,  the  radiating  disks,  the 
pressure  compensation  device,  and  the  radiation  loading.  Although  the  equations  for  this 
transducer  are  more  complicated  than  those  for  the  simple  linear  transducer,  the  form 
of  the  equations  and  the  general  solution  procedure  are  essentially  the  same.  Computed 
results  for  a  single  HLF-6A  transducer  in  water  are  compared  with  computed  results  for 
the  Hydroacoustics,  Inc.  (HAI)  model  when  several  drive  levels  are  used  and  when  the 
input  signal  is  sinusoidal  in  time.  The  agreement  is  quite  good  for  all  drive  levels,  but 
there  are  some  small  departures  at  the  higher  drive  levels  (less  than  1  dB),  which  we 
attribute  to  differences  in  the  radiation  models.  Results  are  also  presented  for  a  proposed 
six-element  array.  Since  the  interelement  spacing  for  this  array  is  quite  large,  the  perfor¬ 
mance  of  each  element  is  not  greatly  different  from  that  of  a  single  element  in  the  free 
field.  Predicted  results  using  our  time-domain  model  are  also  compared  with  experimen¬ 
tal  residts  for  a  three-element  array  which  is  fairly  closely  packed.  The  comparison  is 
reasonably  good  both  in  level  and  in  trends.  Section  4  is  a  summary  of  our  conclusions 
and  recommendations. 
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2  LINEAR  MODEL 


In  this  section  a  time-domain  radiation  model  will  be  developed  and  applied  to  a  simple 
linear  transducer.  The  resulting  linear  model  will  be  used  to  validate  the  time-domain 
approach  to  the  acoustic  radiation  loading  which  will  later  be  applied  to  the  nonlinear 
HLF-6A  hydroacoustic  transducer.  The  simple  linear  transducer  is  shown  in  figure  1.  It 
consists  of  two  rigid  circular  pistons  of  mass  M  placed  symmetrically  at  the  ends  of  a 
rigid  cylindrical  housing.  A  spring  of  stiffness  K  restrains  each  piston  and  each  piston  is 
driven  by  the  same  prescribed  force  F. 

Subsection  2.1  gives  the  mass/spring  relations  for  this  simple  linear  transducer  as  well 
as  a  time-domain  formulation  for  acoustic  loading  on  the  transducer.  The  result  is  a 
system  of  coupled  differential-integral  equations  with  specified  initial  conditions.  The 
solution  procedure  for  solving  this  system  is  described  in  subsection  2.2.  In  subsection  2.3 
the  resxilts  for  a  three-element  array  are  compared  with  those  from  a  frequency- domain 
model  after  the  initial  transients  have  died  out,  in  order  to  validate  the  time-domain 
formulation. 


Figure  1.  Simple  linear  transducer  element. 


2.1  BASIC  EQUATIONS 


The  differential  equation  describing  this  simple  device  is 

Mu  +  Ku  =  F-  F"*-*  (1) 

where  F’’®^  is  the  acoustic  radiation  force  on  the  piston  and  u  is  the  displacement  of  the 
piston.  The  acoustic  radiation  force  is  given  by 

F^^  =  J^pdA  (2) 
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where  p  is  the  acoustic  pressure  and  A  is  the  surface  of  the  piston  in  contact  with  the 
water.  Now  consider  a  planar  array  of  N  of  these  elements.  The  radiation  force  on  the 
m-th  element  of  an  array  can  be  expressed  as  a  linear  function  of  the  piston  normal 
velocities.  In  the  frequency  domain  this  relation  takes  the  form 

^r'M  =  i;^n(u-KH  (3) 

n=l 

where  is  the  radiation  force  on  the  m-th  element,  is  the  mutual  radiation 
impedance  between  the  m-th  and  n-th  elements,  and  u„  is  the  velocity  of  the  n-th  piston. 
These  acoustic  interaction  impedances  sxe  obtained  by  using  the  boundary  element 
program  CHIEF,  which  is  based  on  the  Helmholtz  integral  equation 

|p(C)  =  +  twpu(a)C?(C,<r)]  dS  (4) 

where  S  is  the  radiating  surface  (including  all  array  elements),  v  is  the  normal  velocity,  p 
is  the  fluid  density,  and  G  is  the  free-space  Grin’s  function.  To  obtain  Zmn  the  integral 
equation  4  is  solved  for  p  when  all  the  piston  velocities  are  set  to  zero  except  for  Vn,  which 
is  set  to  1.  The  resulting  pressure  is  integrated  over  the  m-th  piston,  yielding 

When  the  interaction  relation  in  equation  3  is  transformed  into  the  time  domain  via 
a  Fourier  transform,  it  talces  the  form 

K‘\‘)  =  +  E  /'  Cm.(T)t.„(t  -  T)  dr  (6) 

n=l  •'0 

where  is  the  Fourier  transform  of  Zmn  for  ^  is  ^ire  real,  high-frequency 

asymptotic  limit  of  Zmmy  and  Cmm  is  the  Fourier  transform  of  Zmm  —  ^t  turns  out  that 
for  numerical  reasons  it  is  better  to  work  with  accelerations  than  with  velocities.  The 
conversion  of  the  interaction  equations  3  and  5  from  velocities  to  accelerations  gives  the 
frequency-domain  relation 


N 


‘'(w)  =  mmn(u>)anH 

nssl 

(6) 

and  the  time-domain  relation 

‘*(0  =  2  /!  Mmn(T)On(t  -  t)  dr 
n=l 

(7) 

where  mm„(w)  =  Zmn{idj)lioj  and 

PmnlO  =  /  Cmn(r)dT  TTl  ^  Tl 

Jo 

(8) 

=  1  Cmm(r)c/r-f-r*  m  =  n 

Jo 

(9) 
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Figure  2  is  a  plot  of  the  impulse  interaction  coefficient  (^12  appearing  in  equation  5  versus 
normalized  time.  This  plot  illustrates  the  reason  for  expressing  the  convolution  integrals 
in  terms  of  acceleration  instead  of  velocity.  The  fact  that  the  functions  Cmn  ^  n) 
alternate  in  sign  over  a  very  short  interval  gives  rise  to  subtraction  errors  when  the 
convolution  integrals  axe  evaluated.  Figures  3-5  show  the  mutual  impulse  functions  fiu, 
/ii2,  and  appearing  in  equation  7  versus  normalized  time.  Notice  that  these  functions 
are  significantly  nonzero  only  over  a  small  time  interval  (nominally  the  time  for  sound  to 
travel  across  both  the  interacting  elements).  To  obtain  the  mutual  impulse  functions  /Xmn, 
the  CHIEF  [3]  program  was  used  to  calculate  the  mutual  radiation  impedances  Zmn  over  a 
broad  frequency  range.  The  mutual  impedances  were  divided  by  iu  and  transformed  to  the 
time  domain  by  using  an  FFT  routine.  The  self-impulse  function  Hmm  has  a  discontinuity 
at  f  =  0.  To  eliminate  the  resulting  high-frequency  oscillation  in  the  computed  /imm»  we 
use  the  fact  that  the  even  extension  of  a  time  function  has  a  Fourier  transform  which  is 
twice  the  real  part  of  the  original  transform.  Thus,  fimm  can  be  computed  by  talcing  the 
FFT  of  twice  the  real  part  of 


Figure  2.  Mutual  velocity  impulse  coefficient  for  a  pair  of  elements  [dfs  =  2). 
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Figure  5.  Mutual  acceleration  impulse  response  /ii3  (d/s  =  2). 

Combining  the  acoustic  interaction  equation  7  with  the  element  equation  1  gives  the 
following  coupled  system  of  differential  and  integral  equations: 

r* 

Mam  +  KUm  =  ^mn(T)o„(t  -  t)  (It  (10a) 

n=l 

Vm  =  am  (10b) 

Urn  =  Vm  m  =  (lOc) 

where  Fm  is  the  prescribed  force  on  the  m-th  element.  Notice  that  there  is  a  set  of 
equations  like  10  for  each  element  in  the  array  and  that  these  sets  of  equations  are  cou¬ 
pled  through  the  convolution  integrals  in  10a  that  involve  histories  of  all  the  element 
accelerations.  It  will  be  assumed  that  the  elements  are  not  moving  before  t  =  0  and  thus 

Um(0)  =  t;m(0)  =  am{0)  =  0  m  =  1, . . . ,  TV  (11) 

The  differential-integral  equations  10  and  the  initial  conditions  11  describe  the  problem 
to  be  solved  numerically. 
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2.2  SOLUTION  PROCEDURE 


In  this  section  a  solution  procedure  will  be  described  for  the  coupled  set  of  integral- 
differential  equations  10.  As  with  systems  of  ordinary  differential  equations,  the  solution 
is  advanced  successively  by  time  steps  of  length  S  beginning  at  the  initial  time  t  =  0.  To 
describe  the  method  it  is  only  necessary  to  specify  how  the  solution  is  advanced  from  a 
time  ir=r6Xo  the  next  time  U+i  =  i,  + 

It  will  be  convenient  to  first  write  the  equation  10a  in  a  somewhat  different  form. 
When  t  is  in  the  interval  [tr,tr+i]  the  integrals  in  10a  can  be  split  as  follows: 

I  A*mn(7‘)a„(t  -T)dT  =  /  /imn(r)a„(t  -  t)  dr  +  [  fimniT)an{t  -  t)  dr  (12) 
Jo  Jo  Jt~tr 

An  element  in  the  array  does  not  feel  the  effects  of  its  neighboring  elements  immediately 
since  there  is  a  time  delay  required  for  sound  to  travel  between  the  elements.  Thus, 
=  0  for  t  <  dnin/c,  where  dmin  is  the  minimum  distance  between  any  two  elements 
and  c  is  the  sound  speed.  It  will  be  assumed  in  what  follows  that  the  time  step  S  is  smaller 
than  dmin/c.  With  this  assumption 

[  timn(r)an(t  -  t)  dr  =  0  m  ^  n  (13) 

Jo 

since  dfc.  When  m  —  n  the  first  integral  in  equation  12  can  be  approximated 

using  the  trapezoidal  rule  as  follows: 

Mmm(T)a„(t  -  t)  dr  =  I  [Atmm(0)a,n(t)  +  Mmm(i  “  ir)am(fr)]  (<  “  tr)  (14) 
It  follows  from  equations  12  to  14  that 

N  ft 

/  /^mn(T)an(<  -  t)  dr  =  |  [/Xmm(0)a,n(<)  +  /imm(<  “  ir)am(^r)]  (<  “  U)  + 

13/  A^mn(r)an(<-r)dr  (15) 

Substituting  15  into  equation  10a  and  solving  for  am  yields 

®m(0  =  M  ■\--H  (0)(<  —  f  )  ~  “  tr)flm(^r)(^  “  ^r)  ~ 

13  /  Mmn(T)o„(t  -  r)  dr  (16) 

C!ombining  16  with  the  last  two  equations  of  10,  we  obtain 
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Wm(i)  —  ^  (0)(<  —  t  )  ~  \i^mm{i  ^r)<lm(tr)(^  ^r) 

52  /  Mmn(T)a„(i  -  t)  dr  (17a) 

Um(t)  =  t;m(t)  (17b) 

Notice  that  the  right-hajad  side  of  equation  17  involves  only  the  two  unknown  functions 
Um,  Om,  and  past  values  of  which  are  known  at  this  stage. 

The  systems  of  equations  17  can  be  solved  in  the  same  manner  as  a  system  of  ordinary 
differential  equations.  The  only  modification  necessary  is  to  update  the  accelerations  Om 
at  the  end  of  each  step  by  using  equation  16.  It  should  be  pointed  out  that  a  time  step 
must  be  completed  for  each  element  in  the  array  before  the  next  time  step  is  addressed. 
It  is  necessary  to  store  the  time  histories  of  the  accelerations  for  each  element  in  the  array 
and  to  update  these  histories  following  each  time  step.  Since  the  Hmn^  ^le  essentially 
zero  outside  a  small  time  interval,  only  the  relevant  portions  of  the  time  histories  need  to 
be  saved. 

The  results  to  be  presented  later  in  this  paper  were  obtained  using  a  second-order 
Runge-Kutta  method  known  as  Heun’s  method  [4].  For  a  system  of  ordinary  differential 
equations 

»(<)  =  f{yV),t)  (i8») 

»(0)  =  yo  (18b) 

This  solution  procedure  advances  the  solution  yn  at  time  tn  to  the  solution  Pn+i  ^  ^inie 
f„+i  by  using 

yn+l  =  yn  +  I  [j/n  +  fiVn  +  vJ,  tn+l)]  S  (19) 

The  algorithm  given  by  equation  19  together  with  the  formulas  for  updating  accelerations 
can  be  used  to  solve  the  differentiai-integral  equations  given  in  equation  17. 


2.3  RESULTS 

In  this  section  we  will  present  computed  results  for  a  single  element  and  for  a  1  x  3 
array  as  shown  in  figure  6. 

The  results  will  be  checked  against  those  obtained  from  an  independent  frequency-domain 
analysis.  With  the  dimensionless  variables  defined  in  table  1,  equation  10  can  be  written 
as  follows: 


Figure  6.  Three-element  array  configuration. 


Mdm  +  iftim 

dvm 

di 

dum 

di 


-  A 

/imn(f  )an(t  -  f)  df 
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Table  1.  Normalized  variables. 


s/c 

time 

drive  force 

=  \<si/c) 

O 

displacement 

pcs* 

=  ^v(s<7c) 

velocity 

ps^ 

mass 

=  ^<^{silc) 

acceleration 

stiihiess 

The  dimensionless  parameters  in  equation  20  were  chosen  to  make  the  low-frequency 
radiation  mass  of  the  circular  piston  equal  to  the  mass  of  the  piston  and  to  make  the 
piston  diameter  one-tenth  of  a  wavelength  at  the  resonance  frequency.  This  leads  to  the 
dimensionless  parameters 

1  27r* 


(21) 


The  dimensionless  forcing  functions  Fm  were  taken  to  have  the  form 

/■„(<)  = 

=  («) 

where  A  is  the  acoustic  wavelength.  In  the  numerical  examples  the  force  amplitudes 
\Fm\lp<^s^  were  taken  to  be  1.  The  spacing  d  between  elements  for  the  array  example 
was  taken  to  be  2s. 

Figures  7-10  show  a  comparison  of  the  frequency  response  of  the  outer  and  middle 
elements  of  the  1x3  array  when  both  a  frequency-domain  approach  and  a  time-domain 
approach  are  used.  The  agreement  is  very  good.  Figure  11  shows  the  solution  of  the 
differential-integral  equations  10  for  a  single  element.  The  steady-state  normalized  ve¬ 
locity  amplitude  for  this  computed  solution  and  the  value  obtained  from  a  steady-state 
frequency-domain  analysis  of  this  problem  were  both  15.5.  Figure  12  shows  the  computed 
time  response  for  the  middle  element  in  the  1x3  array.  Due  to  the  large  number  of 
oscillations  in  the  given  time  interval  only  the  envelope  of  the  response  is  discernible.  It 
takes  a  very  long  time  for  the  response  to  reach  steady  state  (about  300  times  as  long  as 
it  takes  sound  to  travel  across  the  array). 

Figures  13-15  show  the  effect  of  adding  an  increasing  amount  of  damping  resistance 
to  the  linear  transducer  model.  The  Q  values  given  are  for  the  element  operating  in  a 
vacuiun  and  are  thus  a  measure  of  the  internal  resistance.  As  the  internal  losses  increase, 
the  level  of  the  response  decreases  and  the  time  to  reach  steady  state  decreases.  The  effect 
of  losses  is  important  since  hydroacoustic  transducers  typically  have  very  high  internal 
losses. 
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gure  11.  Time  response  of  a  single  element  in  the  free  field  (s/A  =  0.110). 


Figure  12.  Middle  element  of  three-element  array,  Q=oo. 
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3  HYDROACOUSTIC  MODEL 


Figure  16  is  an  idealized  representation  of  the  internal  components  of  a  hydroacoustic 
transducer  like  the  HLF-6A.  The  transducer  consists  of  two  metallic  radiating  disks  which 
are  clamped  around  their  circumference  to  a  cylindrical  housing.  The  disks  are  driven 
at  their  centers  by  two  drive  pistons,  which  in  turn  are  driven  by  the  pressure  in  an  oil- 
filled  drive  cavity.  The  drive  cavity  pressure  is  modulated  by  a  sliding  valve  that  controls 
the  flow  into  and  out  of  the  cavity.  This  valve  is  driven  by  a  hydraulic  servovalve.  The 
servovalve  is  called  the  first-stage  valve  and  the  main  valve  is  called  the  second-stage  valve. 
As  the  second-stage  valve  slides  back  and  forth,  the  drive  cavity  is  alternately  coupled  to 
the  supply  and  return  ports  of  a  hydraulic  pump.  The  pump  is  designed  to  maintain  a 
constant  pressure  at  the  supply  port.  The  supply  pressure  is  much  higher  than  the  return 
pressure. 

This  section  contains  descriptions  of  the  mathematical  models  for  the  main  compo¬ 
nents  of  the  HLF-6A  transducer.  The  equations  for  the  first-stage  servovalve  are  given  in 
subsection  3.1  and  the  equations  for  the  second-stage  valve  are  given  in  subsection  3.2. 
The  equations  for  the  added  stiffness  due  to  the  pressure  compensation  mechanism  are 
given  in  subsection  3.3.  The  mathematical  model  for  the  raxliating  disks,  drive  cavity,  and 
the  acoustic  loading  is  presented  in  subsection  3.4.  The  model  for  the  acoustic  loading  is 
a  straightforward  extension  of  the  model  previously  described  for  the  simple  linear  trans¬ 
ducer.  A  sximmary  of  the  differential-integral  equations  for  an  array  of  hydroacoustic 
elements  is  presented  in  subsection  3.5.  Subsection  3.6  presents  results  obtained  from 
this  model  for  both  single  elements  and  arrays.  Single-element  results  are  compared  with 
those  of  an  existing  HAI  model  for  various  drive  levels.  Results  for  a  proposed  six-element 
array  are  compared  with  those  of  a  single  element.  Model  results  for  a  close-packed  three- 
element  array  are  compeired  with  experimental  results. 


3.1  FIRST-STAGE  SERVOVALVE 


This  section  describes  a  mathematical  model  for  the  first-stage  servovalve  used  in  the 
HLF-6A  and  its  associated  drive  circuit.  A  diagram  of  the  servovalve  and  drive  circuit  is 
shown  in  figure  17. 

The  internal  components  of  the  servovalve  (enclosed  by  a  dotted  box  in  figure  17)  were 
not  modeled  in  detail.  Instead  an  empirical  relationship  of  the  form 

(23) 

was  used  to  relate  the  drive  current  /  and  the  servovalve  displacement  x«,.  We  will  next 
show  how  to  relate  the  servovalve  displacement  to  the  servovalve  output  flow  Q„  and  the 


Figure  16.  Idealized  hydroacoustic  transducer. 
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Figure  17.  First-stage  servovalve  and  associated  drive  circuit. 


(24) 


main  valve  displacement  i„.  Let  us  define  Pi  and  as  follows: 

p  _[P,  if  >  0  p  _  /  -Pr  if  Xw  >  0 
^~\Pr  ifi«,<0  "  IP,  ifx.„<0 

Then  the  flow  Qt,  through  the  servovalve  satisfies  the  equations 

Pi-p;  =  p^-P2  =  i,^  (25) 

®;t» 

where  is  the  sign  of  and  is  a  constant  that  depends  on  the  geometry  of  the  valve 
and  the  density  of  the  fluid.  Thus 

IvQv  >  0  (26) 

The  flow  Qv  also  satisfies  the  orifice  equations 

pi-p.=p',-pi= (27) 


where  ai  depends  on  the  geometry  of  the  two  identical  output  orifices  and  on  the  density 
of  the  fluid.  The  flow  Qp  through  the  shunt  orifice  satisfies 

P,-P:  =  I,a,Ql  (28) 

where  /,  is  the  sign  of  Qp  and  aj  is  a  constant  that  depends  on  the  geometry  of  the  shunt 
orifice  and  on  the  density  of  the  fluid.  The  flows  Q*,  and  Qp  also  satisfy  the  balance 
equation 

Q.  =  Q,  +  Q.  (29) 

It  follows  from  equation  29  that 

l,Qp  =  l,iQv-Qe)>0  (30) 

Addition  of  the  equations  25,  27,  and  28  yields  the  equation 

Pi-P2  =  I,ot2Ql  +  2haiQl  +  24^  (31) 

Combining  equations  29  and  31  gives 

Pi  -  P2  =  iMQv  -  Qcf  +  2hotiQl  +  2/v  ^  (32) 

or 

^/,a2  +  2hai  +  Ql  -  i2I,a2Qe)Qv  +  (/.asQ^  -  Pi  +  Pj)  =  0  (33) 

The  quadratic  equation  33  has  the  solutions 


-B  +  Iry/B^  -  4 AC 
2A 


(34) 
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where  Ir  —  d:l  and 


A  —  I,  +  2/v  fai  + 

(35a) 

B  =  —21,Q2Qe 

(35b) 

C  =  J,a,Ql-Pt  +  P2 

(35c) 

Since  must  be  real,  we  must  have 

-  iAC  >  0 

(36) 

As  we  are  ignoring  the  compliance  of  the  end  cavities,  the  flow 
displacement  x^  by 

is  related  to  the 

Qe  ~  A^X„ 

(37) 

FVom  equations  28,  29,  and  37  it  follows  that  the  driving  force  on  the  second-stage 
valve  spool  is  given  by 


Fv  =  (Pe  -  Pe)Av  =  -  Ar,Xy)  (38) 

However,  can  be  determined  by  solving  the  quadratic  equation  33,  given  x„  and  hence 
Qe*  The  quantities  /«,  and  Ir  are  determined  so  that  the  conditions  26,  30,  and  36  are 
satisfied.  Thus,  we  have  shown  that 


F„=f{x„,x„,t)  (39) 

3.2  SECOND-STAGE  VALVE 


The  second-stage  valve  is  pictured  in  figure  18.  The  valve  spool  is  shown  in  the  neutral 
position  {xv,  =  0).  The  overlap  between  the  valve  spool  and  the  valve  sleeve  on  each  side 
of  the  drive  cavity  opening  when  the  valve  spool  is  in  the  neutral  position  is  called  the 
valve  overlap  (wol). 


The  equations  describing  the  flow  through  the  supply  (return)  port  of  the  valve  into 
the  drive  cavity  axe  different  when  the  supply  (return)  port  is  open  and  when  it  is  closed. 
When  the  supply  (return)  port  is  open,  the  flow  at  low  frequencies  is  given  approximately 
by  the  flow  through  an  orifice  whose  area  is  changing  slowly  with  time.  It  is  shown  in 
a  previous  paper  [1, Appendix  D]  that  the  pressure  drop  across  the  open  supply  port  is 
given  by 


n  n  ^  Ph  Qt\Qi\ 


(40) 


where  Pj  is  the  pressure  on  the  upstream  side  of  the  supply  port,  Pj  is  the  pressure  on 
the  drive  cavity  side  of  the  supply  port,  is  the  area  of  the  supply  port  opening,  and 
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Figure  18.  Second-stage  valve. 

7  is  the  ratio  of  the  flow  stream  area  at  the  point  of  greatest  contraction  to  the  valve 
opening  area.  The  pressure  Pi  is  different  from  the  supply  pressure  P,  since  there  are 
passages  and  fittings  between  the  pump  and  the  valve.  Similarly,  the  pressure  drop  across 
the  return  port  when  it  is  open  is  given  by 


P3-P3  = 


Ph  Qt\Qt\ 


(41) 


2 

where  P3  is  the  pressure  on  the  downstream  side  of  the  return  port  and  Arp  is  the  area  of 
the  return  port  opening. 


We  next  derive  expressions  for  the  areas  A,p,  Arp  as  a  function  of  valve  position.  The 
supply  port  opening  is  pictured  in  figure  19.  The  rounded  corners  shown  on  the  valve 
sleeve  and  the  valve  spool  are  greatly  exaggerated  (rad=0.002  in.).  When  the  supply  port 
is  open  we  have 

x„  =  \l {xv$  —  wol  -I-  2  •  rad*  -|-  (rcl  -f-  2  •  rad)^  —  2  •  rad  (42) 

where  is  the  displacement  of  the  valve  from  neutral  position  (a  positive  increase  in 
tends  to  open  the  supply  port).  The  area  A»p  of  the  supply  port  opening  is  given 
E^proximately  by 

Atp  =  irDyXy  (rcl  <  Zl„)  (43) 

where  is  the  diameter  of  the  valve  spool. 

Similarly,  when  the  return  port  is  open  we  have 

Xv  —  \/  (— Xv*  —  wol  -H  2  •  rad*  +  (rcl  -I-  2  •  rad)*  —  2  •  rad  (44) 

and 

Arp  =  7rZ?„Xv  (rcl  <  Dv).  (45) 
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Figure  19.  Supply  port  geometry. 


(46) 

(47) 


The  pressure  drops  P,  —  Pi  and  P2  —  Pp  are  ©ven  by 

D  D  _  r  /S  I  PhQ»\Qt\ 


P$  —  Pi  =  LfQs  + 


27MJ 


P2  —  Pd  =  L,oQ»  + 


PhQ,\Q$\ 


*  ^  '  2ri^Al„  "  " 

where  L,  is  the  supply  flow  path  inertance,  //«,  is  the  supply  output  flow  path  inertance,  A, 
is  the  supply  flow  path  equivalent  area,  and  A,o  is  the  supply  output  flow  path  equivalent 
area.  The  pressure  drops  Pp  —  P2  and  P3  —  Pr  axe  given  by 


P3  —  Pr  —  LtQt  + 


^^QrlQrl 

2»7M3 


and 

Pd  -  ft  =  L,.Q.  +  ^1^^  (61) 

where  Lr  is  the  return  flow  path  inertance,  Lro  is  the  return  output  flow  path  inertance,  Ar 
is  the  return  flow  path  equivalent  area,  and  Aro  is  the  return  output  flow  path  equivalent 
area. 
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Combining  equations  40,  48,  and  49  yields 


+  P.  Pd  2^1 


(52) 


when  the  supply  port  is  open.  Similarly 

(i.  +  -  ft  -  0  (i  +  ^  +  -5^)  9,W,|  (53) 

when  the  return  port  is  open. 

When  the  supply  port  is  closed  the  leakage  flow  is  given  by  [l,p.l5]: 

Q,  =  ,gii(ft,  -  Pd)  [- +  V^t2«i  +  t3|P.-PDl]  (54) 

where 

o  12fiXL  ,2  2i7’ir*£)2(rcl)*  , 

= ;f5;:(;d)f  - Tk — ■  (“) 

and  /X  is  the  shear  viscosity  coeflKcient  of  the  fluid.  Similarly,  when  the  return  port  is 
closed 

Qr  =  Sgn(Px,  -  Pr )  [- \1^Rl  +  >J\^RL^l^\PD-Pr\ ]  (56) 

where  now  xx,  =  x„,  +  wol. 


The  equation  of  motion  for  the  second-stage  valve  spool  can  be  written  in  the  form 

■Ww«x„,  -}-  2KfXv»  =  Pv  +  Pjiow  (57) 

or  eqmvalently 

Vv$  ~  ~77  (.Pv  "b  P flow  ~~  2P»Xv»)  (58a) 

x„,  =  Vvt  (58b) 

where  Xv«  is  the  lineal  valve  displacement,  Mv«  is  the  mass  of  the  valve  spool,  K,  is  the 
stiffness  of  each  of  the  two  restraining  springs,  Pv  Is  the  driving  force  on  the  spool  given 
in  equation  38,  and  Pjiow  is  the  axial  force  on  the  spool  due  to  the  flow  through  the  valve. 
The  force  P/iou,  is  given  approximately  by  [5]: 

Pjiow  =  -8sa{x„t)  pkQU  cos  9  (59) 

where  pk  is  the  density  of  the  hydraulic  fluid,  Q  is  the  flow  through  the  valve  (supply  or 
return),  U  is  the  velocity  of  the  jet  at  the  point  of  maximum  contraction  (vena  contracta)^ 
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and  CO8  0  =  0.36  ($  is  the  angle  of  the  flow  relative  to  the  axis).  Since  Q  —  rjU A,  where 
A  is  the  area  of  the  valve  opening,  equation  59  can  be  written 

Fjiou,  =  -0.36sgn(x„,)^^  (60) 

The  area  of  the  valve  opening  is  given  approximately  by 

A  =  tDvXv  (61) 

Thus,  equation  60  can  be  written 

= -0.36sgn(x„,)~^—  (62) 

Since  the  flow  through  the  closed  port  is  small,  we  can  write  approximately 

F,i„  =  -0.36  sgn(*..)  (63) 

3.3  PRESSURE  COMPENSATION 

In  order  to  allow  the  transducer  to  operate  at  large  depths  the  interior  of  the  transducer 
is  filled  with  a  pressurized  gas.  The  pressure  of  the  gas  is  adjusted  to  be  equal  to  the 
hydrostatic  pressure  in  the  water  at  the  operating  depth.  Thus,  the  pressure  Pq  of  the 
interior  gas  is  adjusted  so  that 

Po  ~  P »tm  +  pgD  (64) 

where  Patm  is  atmospheric  pressure,  p  is  the  density  of  the  water,  g  is  the  acceleration  of 
gravity,  and  D  is  the  operating  depth. 

In  addition  to  providing  depth  compensation,  the  pressurized  gas  also  presents  a  dy¬ 
namic  stiffness  to  each  of  the  radiating  disks.  By  symmetry  it  is  only  necessary  to  consider 
one  of  the  radiating  disks  and  half  of  the  interior  volume.  When  the  disk  vibrates  the 
interior  volume  changes  which  in  turn  produces  a  change  in  pressure.  It  will  be  assumed 
that  the  pressurized  gas  obeys  the  ideal  gas  law 

PV'f  =  PoVS'  7  =  1-4  (65) 

where  7  is  the  adiabatic  gas  constant,  P  is  the  pressure,  V  is  the  volume  of  gas,  Po  is  the 
initial  pressure  charge,  and  Vo  is  the  initial  volume  of  gas.  Let  P  and  V  be  the  changes 
in  pressure  and  volume  produced  by  the  vibration  of  the  disk,  i.e. 

P  =  Po  +  P  V  ^  Vo +  V  (66) 
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Combining  equations  65  and  66,  we  obtain 

Assuming  that  P  ^  Pq  and  V  Vo,  it  follows  from  equation  67  that 

or 

P  =  -7^^  (69) 

If  u  is  the  normal  displacement  of  the  disk,  then  the  volume  change  V  is  given  by 

V^j^udA  (70) 

The  displacement  u  can  be  approximated  by 


U  =  Ui^i  -f  +  **4^4 

(71) 

where  and  ^4 

unity  at  the  disk  center). 

are  the  dominant  mode  shapes  of  the  disk  (normalized  to  be 
Combining  equations  69-71,  we  obtain 

‘'0  n=l 

”  y  1/  OlnUn 

^0  nsl 

(72) 

where  On  is  defined  by 

II 

>>^1  *- 

(73) 

The  finite-element  modal  force  component  Fm  is  given  by 

Fm  =  J  =  FAOm 

(74) 

Thus 

p  _ 

^ m  —  ^  OCfnOlnXiyi 

^  Ksl 

nsl 

(75) 

where 

p  -  y^oA\ 

'^mn  j, 

Vo 

(76) 
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Let  fi  be  the  fraction  of  the  interior  volume  that  is  filled  with  gas.  Then 


Vo  =  0Ah  (77) 

where  A  is  the  area  of  the  disk  and  h  is  the  half  height  of  the  cylindrical  housing.  Sub¬ 
stituting  equation  77  into  equation  76,  we  get 

(78) 

ph 

Equation  75  is  the  desired  dynamic  stiffness  relation. 


3.4  RADIATING  DISKS  AND  DRIVE  CAVITY 


The  radiating  disks  were  modeled  by  using  the  finite-element  method.  Pour  symmetric 
modes  were  used  to  represent  the  motion  of  the  disks.  The  modes  were  normalized  so  that 
the  normal  component  of  displacement  at  the  disk  center  was  unity.  In  terms  of  these 
modes  the  equation  of  motion  for  the  m-th  element  can  be  written 

+  =  •  =  1,2,3, 4  (79) 

where  Mi  is  the  modal  mass  of  the  t-th  mode,  Ki  is  the  modal  stiffness  of  the  t-th  mode, 
is  the  t-th  modal  displacement  of  the  m-th  element.  Pm  is  the  drive  cavity  pressure  of 
the  m-th  element,  Ap  is  the  area  of  the  drive  piston,  and  is  the  t-th  modal  component 
of  the  acoustic  radiation  force  on  the  m-th  element.  This  equation  can  be  written  in  the 
alternative  form 


MicUt)  +  KiuUt)  =  P„(4).4,  -  F;(()  (80a) 

<(.*)  =  <C(‘)  (80b) 

»=(<)  =  <(<)  i=l,2,3,4  (80c) 

Equation  7  for  the  radiation  force  in  the  linear  model  can  be  extended  to  this  problem  as 
follows: 

=  1  =  1, 2, 3, 4  (81) 

j=l  nsl 

Combining  equations  80  and  81,  we  obtain 

Mi0.l{t)  +  K>uUt)  =  (82a) 

jsl  jsl  r»arl 

iUt)  =  «i.(<)  (82b) 

0i.(*)  =  <(«)  i  =  l,2,3,4  (82c) 
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In  order  to  solve  the  differential-integral  equations  for  this  problem,  we  will  rewrite 
equations  82  as  was  done  in  the  linear  problem.  Let  t  be  in  the  interval  U  <t  <  tr^x. 
Then  the  integrals  on  the  right-hand  side  of  equation  82  can  be  split  as  follows: 

/  -  r)  dr  =  /  '  -  t)  dr  -f  /  -  t)  dr  (83) 

JO  Jo  Jt-tr 

Since  an  element  does  not  feel  the  effects  of  its  neighboring  elements  immediately,  it 
follows  that 

M‘mni‘r)a'i{t-T)dT  ==0  m^n  (84) 

When  m  =  n  the  first  integral  in  equation  83  can  be  approximated  by  using  the  trapezoidal 
rule  as  follows: 

-  r) dr  =  (/C.(0K W  +  -  <.)<(<.)]  (8S) 

Inserting  equations  83-85  into  equation  82  yields 

^  ^  j  jsX  jsl 

jsl  nsl  J*-tr 

(^■^)  2  /“iLC*  -  ‘r)aj.(‘.)  (86») 

<(t)  =  <(t)  (86b) 

“UO  =  oUO  i  =  l, 2,3.4  (86e) 

The  three  equations  in  86  can  be  solved  for  the  a)^  in  terms  of  the  remaining  variables. 


The  pressure  Pm  in  the  drive  cavity  satisfies  the  equation 

CPm(t)  =  net  flow  into  the  cavity  =  Q^(t)  -  ApV^{t)  (87) 

where  C  is  the  compliance  of  the  cavity,  is  the  center  velocity  of  element  m,  and 
is  the  flow  into  the  drive  cavity  of  the  m-th  element  through  the  second-stage  valve.  Since 
the  modes  were  normalized  to  unity  at  the  center,  the  velocity  can  be  written 

«r(‘) = !:»:.(«)  (88) 

i=\ 

The  net  flow  is  given  by 

=  QUt)  -  QmCO  (89) 


where  and  are  the  supply  and  return  flows  for  the  m-th  element.  Substituting 
equations  88  and  89  into  equation  87  produces 


«»(0  =  5 


t=l 


(90) 
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3.5  SUMMARY  OF  EQUATIONS 


In  this  subsection  we  will  present  a  summary  of  the  differential-integral  equations  for  a 
hydroacoustic  array  obtained  in  previous  subsections.  The  equations  will  be  given  for  the 
m-th  element  of  an  array.  There  will  be  a  system  of  equations  like  this  for  each  element 
in  the  array. 

Radiating  disk  equations: 


u*  =  a* 

fft  ffi 


“m  = 


*  =  1,2,3, 4 


where 


+  Kiuln  =  PmAp  -  ^  kijui,  f  “  r)  dr 


Drive  cavity  equation: 


t=l 


Valve  spool  equations: 


where 


+  Ffiov,  —  2K,Xv,) 

iVlv< 

Uti# 


Fv  — 

Ffiou.  =  -0.36sgn(x,,)^'''-^l'^^*^ 

T/7rZ/„i„ 


Valve  flow  equations: 


(Lt  +  Lto)Q$  =  P$~  Pd 


=  P.  -  Pn  -  ^  f  -L 


2*7* \Ai 


ir^Dlxl 


) 


Qt\Qi\  supply  port  open 
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{Lr  +  Lro)Qr  =  "  Pr  "  ^  ^  <2r|Qrl  rctuni  port  Open 


Q,  =  8gn(P.  -  Pd)  [- +  \I\^Rl-^^\P,- Pd\^  supply  port  leaJcage 

Qr  =  8gn(PD  -  Pr)  [“  +  ^|Pd  -  P- 1 ]  retum  port  leakage 

The  initial  conditions  for  u)„  and  (t  =  1,2, 3, 4),  Xm*  u^,  and  are  taken  to 
be  zero.  The  initial  condition  on  the  drive  cavity  pressure  Pm  is  taken  to  be  (P«  +  Pr)!^- 
The  solution  procedure  as  outlined  for  the  linear  problem  is  applicable  to  the  solution  of 
this  set  of  differential-integral  equations. 


3.6  RESULTS 


This  subsection  contains  computation  results  from  the  time-domiun  model  for  both 
single  elements  and  arrays  of  HLF-6A  transducers.  The  following  normalizations  have 
been  employed  in  the  figures: 


1.  FVequency  has  been  shifted  by  a  value  designated  as  FL. 

2.  Time  has  been  divided  by  s/c,  where  s  is  the  diameter  of  an  element  and  c  is  the 
soxmd  speed. 

3.  All  so\irce  levels  have  had  the  single-element  target  source  level  subtracted. 

4.  Array  source  levels  also  have  had  20  logjo(A^)  subtracted,  where  N  is  the  number  of 
elements. 

5.  Velocities  have  been  divided  by  the  sound  speed  c. 

6.  Drive  cavity  pressures  have  been  divided  by  the  supply  pressure. 


Figures  20-23  show  a  comparison  of  results  from  the  NRaD  model  and  the  HAI  model 
for  sinusoidal  drives  of  four  different  levels  and  a  number  of  different  frequencies.  At  each 
drive  frequency  the  fundamental  component  of  the  time  response  is  plotted.  A  drive  levd 
of  0  dB  corresponds  to  half  the  maximum  displacement  of  the  servovalve.  For  all  drive 
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levels  the  agreement  is  quite  good,  but  there  are  small  differences  (less  than  1  dB)  at  the 
higher  drive  levels.  HAI  models  the  load  (disk+pressure  compensation+radiation)  by  a 
mass,  stiffness,  resistance,  and  shape  factor  for  each  drive  frequency.  The  values  of  these 
parameters  axe  chosen  to  match  the  drive  point  impedance  at  the  drive  frequency.  This 
load  is  not  correct  at  harmonics  of  the  drive  frequency  produced  by  the  nonlinearities 
of  the  device.  In  addition  this  approach  is  not  applicable  to  nonsinusoidal  drive  signals. 
Figures  24-27  show  the  same  comparison  as  above  except  that  the  NRaD  model  was 
modified  to  use  the  HAI  loads.  For  all  drive  levels  the  results  are  virtually  identical.  This 
leads  us  to  believe  that  the  differences  in  figures  20-23  axe  due  to  differences  in  the  way 
the  load  was  modeled. 

The  next  eight  figures  contain  comparisons  related  to  the  drive  point  impedance  of 
the  load.  Figures  28  and  29  show  a  comparison  of  the  magnitude  and  phase  of  the  drive 
point  impedance  between  the  NRaD  time-domain  model  and  the  NRaD  frequency-domain 
model.  In  the  time-domain  model  only  the  fundamental  components  were  used  to  form  the 
impedance.  The  agreement  is  quite  good  over  the  entire  frequency  band.  Figures  30  and 
31  show  the  same  comparison  for  the  HAI  time-domain  and  frequency-domain  models. 
Again  the  agreement  is  excellent.  Figures  32  and  33  show  a  comparison  of  the  drive  point 
impedance  magnitude  and  phase  between  the  NRaD  frequency-domain  model  and  the 
HAI  frequency-domain  model.  The  agreement  is  excellent.  Figures  34  and  35  show  the 
same  comparison  for  the  time-domain  models.  Again  the  agreement  is  good.  Since  the 
NRaD  and  HAI  models  produce  virtually  the  same  drive  point  impedance  at  the  load, 
it  appears  that  any  differences  in  results  must  be  due  to  higher  harmonics  of  the  drive 
frequency. 

The  next  four  figures  show  results  for  a  three-element  line  array  in  which  the  center-to- 
center  spacing  is  twice  the  element  diameter.  Figure  36  shows  a  comparison  of  normalized 
source  level  computed  with  the  NRaD  time-domain  model  and  the  experimental  results 
obtained  at  Lake  Seneca  by  Naval  Research  Laboratory/Underwater  Sound  Reference 
Detachment  (NRL/USRD).^  The  comparison  is  reasonably  good  both  in  levels  and  in 
trends,  but  there  appears  to  be  a  small  frequency  shift  in  the  peaks.  Figure  37  shows  a 
comparison  of  model  results  for  a  single  element  and  for  the  three-element  array.  This 
comparison  indicates  that  at  this  element  spacing  the  array  source  level  can  not  be  ob¬ 
tained  by  adding  20  logiQ(3)  to  the  single-element  source  level.  Figures  38  and  39  show 
the  time  response  of  the  center  velocity  and  drive  cavity  pressure  for  the  middle  element 
of  the  three-element  array.  The  fact  that  the  transient  behavior  dies  out  very  quickly  is 
due  to  inherently  high  losses  in  the  second-stage  valve  flow. 

The  next  three  figures  are  for  a  proposed  six-element  line  array  in  which  the  center-to- 
center  spacing  is  9.26  times  the  element  diameter.  Figures  40  and  41  show  a  comparison 
of  the  time  response  of  the  center  velocity  and  drive  cavity  pressure  for  a  single  element 

^Data  supplied  by  Steven  Black  of  NRL/USRD  in  a  private  communication 
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versus  a  middle  element  of  the  six-element  array.  At  this  large  element  spacing  the  results 
are  virtually  identical.  Figure  42  shows  that  the  source  level  for  this  six-element  array 
can  be  obtained  by  adding  20  logiQ(6)  to  the  single-element  source  level. 
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Figure  28.  NRaD  impedance  magnitude  comparison. 
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Figure  29.  NRaD  impedance  phase  comparison. 
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Figure  30.  HAI  impedance  magnitude  comparison. 
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Figure  31.  HAI  impedance  phase  comparison. 
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Figure  40.  Velocity  comparison,  middle  element  of  1  X  6  array  vs  single  element. 
U3  ...n-i  SINGLE  ARRAY 
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Figure  41.  Drive  cavity  pressure  comparison,  middle  element  of  1  x  6  array  vs  single 
element. 


2 
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Figure  42.  Normalized  source  level,  single  element  vs  six-element  array. 


4  CONCLUSIONS  AND  RECOMMENDATIONS 


A  mathematical  model  has  been  derived  for  an  array  of  hydroacoustic  transducers 
and  applied  to  the  HLF-6A.  The  major  new  feature  of  this  model  is  the  time-domain 
approach  to  the  array  acoustic  interaction  problem.  This  approach  has  been  validated  for 
an  array  of  simple  linear  transducers  by  comparison  with  a  frequency-domain  approach. 
This  linear  problem  was  also  used  to  study  the  effects  of  internal  damping  on  the  time 
required  to  reach  steady  state.  It  was  found  that  when  the  damping  was  small  the  time 
to  reach  steady  state  could  be  very  long  compared  with  that  of  a  single  element.  As 
the  damping  was  increased,  the  time  to  reach  steady  state  decreased  dramatically  and 
the  steady  state  level  became  lower.  Results  for  a  single  HLF-6A  transducer  have  been 
compared  with  corresponding  results  from  HAPs  computer  model.  The  agreement  is  in 
general  very  good.  Slight  differences  are  attributable  to  differences  in  the  way  acoustic 
radiation  loading  is  handled.  The  results  of  our  time-domain  model  were  also  compared 
with  experimental  results  for  a  close-packed  three-element  array.  The  comparison  showed 
the  same  general  trends  and  levels,  but  appeared  to  be  slightly  shifted  in  frequency.  A 
proposed,  wide-spaced,  six-element  array  was  also  modeled.  At  this  proposed  spacing  the 
interaction  effects  were  minimal.  In  both  the  three-  and  six-element  arrays  the  time  to 
reach  steady  state  was  about  the  same  as  for  a  single  element  since  the  internal  damping 
of  the  valve  is  large. 

In  the  future  we  would  like  to  investigate  an  alternative  approach  in  which  the  in¬ 
teractions  are  modeled  at  the  drive  point  impedance  level  instead  of  at  the  radiation 
impedance  level.  This  would  avoid  the  difficulty  of  choosing  which  modes  to  incorporate 
in  the  model.  We  would  also  like  to  further  modularize  our  computer  code  so  that  other 
differential  equation  solvers  could  be  employed. 
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